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Abstract. We show that the laws of Zipf and Benford, obeyed by scores 
of numerical data generated by many and diverse kinds of natural phe- 
nomena and human activity are related to the focal expression of a gen- 
eralized thermodynamic structure. This structure is obtained from a de- 
formed type of statistical mechanics that arises when configurational 
phase space is incompletely visited in a severe way. Specifically, the re- 
striction is that the accessible fraction of this space has fractal properties. 
The focal expression is an (incomplete) Legendre transform between two 
entropy (or Massieu) potentials that when particularized to first digits 
leads to a previously existing generalization of Benford's law. The inverse 
functional of this expression leads to Zipf 's law; but it naturally includes 
the bends or tails observed in real data for small and large rank. Re- 
markably, we find that the entire problem is analogous to the transition 
to chaos via intermittency exhibited by low-dimensional nonlinear maps. 
Our results also explain the generic form of the degree distribution of 
scale-free networks. 
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1 Introduction 

Over more than half a century, observers of the astonishing ubiquity of the 
empirical laws of Zipf and Benford have been puzzled by their seeming universal 
validity and intrigued about the plausible answer to the central question of why 
they appear in so many contexts. As it is widely known, Zipf's law refers to the 
(approximate) power law that is displayed by sets of data (populations of cities, 
words in texts, impact factors of scientific journals, etc.) when these are given 
a ranking (in relation to size of populations, frequency of words, magnitude of 
impact factors, etc.) [1]. Benford's law is a well-known simple logarithmic rule 
for the frequency of first digits found in listings of data (stock market prices, 
census data, heat capacities of chemicals, etc.) [2]. 

It has been argued that Benford's law is a special case of Zipf's law [3]. Indeed 
the relationship between the two has been made explicit some years ago [4] by 
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first obtaining a generalization of Benford's law from the basic assumption that 
the underlying probability distribution P(N) of the data N under consideration 
is scale invariant and therefore has the form of the power law P(N) <~ N~ a , 
a > 1. A simple integration over P(N), to obtain the relative probability for 
consecutive integers n and n + 1, leads, when a = 1, to 7r(n) = log(l + 
which is Benford's law. The next step in Ref. [4] was to obtain the rank k from 
P(N), this time as an integration over P(N) from N(k), the number of data 
that define the rank k, to a finite number iV max that corresponds to the first 
value of the rank k. In the limit iV max — > oo they obtain N(k) ~ fcVC 1 -") that 
is Zipf's law with exponent 1/(1 — a) when a > 1. For many sets of real data 
a ~ 2 and the standard Zipf law is a = 2 [4]. 

Here we expand on the results of Ref. [4]. Our first, simple, step is to keep 
JV max finite, but as we argue below, this consideration facilitates the articulation 
of a major inference on the physical nature of the laws of Zipf and Benford. 
We contend that these laws are related to a general thermodynamic expression, 
albeit for a special type of thermodynamic structure obtained from the usual 
via a scalar deformation parameter represented by the power a. The general 
thermodynamic expression is seen to represent an (incomplete) Legendre trans- 
form (akin to a Landau free energy or a free energy density functional) between 
two thermodynamic potentials, the expression relating the corresponding parti- 
tion functions becomes a generalized Zipf's law. We identify these quantities as 
well as the conjugate variables involved, which are the rank k and the inverse 
of the total number of data A/ --1 . We also reason that this kind of deformed 
thermodynamics arises from the existence of a strong impediment in accessing 
configurational phase space, that materializes in only a fractal or multifractal 
subset of this space being available to the system. A quantitative consequence 
of considering iV max finite is the reproduction of the small-rank bend displayed 
by real data before the power-law behavior sets in. The power law regime in the 
theoretical expression persists up to infinite rank k — > oo, indicating a kind of 
'thermodynamic limit'. We illustrate this feature by comparison with on hand 
data for frequencies of English words in texts [5] . We refer to the application of 
this scheme to the degree distribution of scale-free networks. 

A subsequent development is the identification of a strict analogy between the 
aforementioned thermodynamic expression and that for (all) the trajectories at 
the transition to chaos via intermittency in nonlinear low-dimensional maps, the 
so-called tangent bifurcation [6]. These critical trajectories follow [7] the exact 
closed form of the functional-composition Rcnormalization Group (RG) fixed- 
point map [8] [6]. Consequently, we associate the same statistical-mechanical 
structure to the nonlinear dynamics at this transition. Further, we examine the 
modifications brought upon the generalized law of Zipf by the corresponding 
shift of the map out of tangency into the chaotic regime. These consist in the 
introduction of an upper bound for the rank k and the reproduction of the tail 
observed for large rank in real data. We illustrate our scheme by comparison 
with available numerical data for the so-called eigenfactor of physics journals 
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[9], industrial production rates [10], and carbon emissions [11]. The analogy 
indicates that the most common value for the index a should be a = 2. 

Lastly, we make use of the statistical-mechanical interpretation to extend our 
analysis. We presuppose that the Legendre transformation expressed by these 
laws can be finalized in the usual way and eliminate the variable Af^ 1 in favor 
of k. To accomplish this step it is necessary to specify the (partition) function 
-Wmax(-A/" -1 ), a feature of the available data or a prerogative of the data collector, 
and evaluate the 'equation of state' k(Af ). In doing this it becomes clear that 
the universality of the laws is due to the general form of the incomplete Legendre 
transformation, while the initial and transformed thermodynamic potentials are 
particular of the data in hand. 

Thus, in the next Section 2 we reproduce the expressions in Ref. [4] relevant 
to our purposes. In Section 3 we describe the generalized statistical-mechanical 
structure we observe in these expressions. In the following Section 4 we present 
the parallelism between the ranking of data and the critical dynamics at the 
tangent bifurcation in nonlinear maps and describe the finite size effect of the 
former. In Section 5 wc extend the statistical-mechanical description and draw 
conclusions on the apparent universality of the aforementioned empirical laws. 
We conclude in Section 6 with a summary and discussion. A partial preliminary 
account of the contents of this paper appeared in Ref. [12] 

2 Derivation of the Laws of Benford and Zipf 

Denote by P(N) the probability distribution associated to the set of data under 
consideration (e.g., the distribution obtained from a histogram generated by 
data - a total of N numbers - giving the magnitudes of the population of a set 
of countries). Under the assumption of scale invariance the distribution has the 
form of a power law P(N) <~ N~ a , a > 0. The probability of observation of the 
first digit n of the number N is given by [4] 

71+1 

7r(n)= [ N- a dN=— i —\(n + l) 1 - a -n 1 - a }, (1) 
J 1 - a 

n 

a^l, from which one obtains Benford's law n(n) = log(l + n^ 1 ) when a = 1. 

The set of Af factual data numbers can be ranked and compared with ranking 
of another set of also Af numbers extracted from the basic distribution P(N) <~ 
N~ a . The rank k is given by k = Af II(N(k) , iV max ) where [4] 

JV m „ 

n(N(k),N max )= I N~ a dN 

N(k) 

= I ±-[N^-N{k) 1 - a ], (2) 

a ^ 1, where N max and N(k) correspond, respectively, to rank k = 0, and 
nonspecific rank k > 0. Solving the above for N(k) in the limit iV max ^> 1 yields 
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Zipf's law N(k) ~ fc 1 ^ 1- "). Eq. (2) introduces a continuum-space variable for 
the rank k in which the first value of the rank is k = 0. This is a departure 
from the usual representation with first rank k = 1 and the following ranks 
given by successive natural numbers. This approach corresponds to a continuum 
variable description suitable for large data sets, and for which restriction to 
integer values of the rank can be obtained by use of suitable values for the lower 
limits of integration N(k) in Eq. (2). 



3 Generalized laws of Benford and Zipf as thermodynamic 
relations 

Consider the g-dcformcd logarithmic function log 9 (a;) = (1 — (j) -1 [a: 1-9 — f] 
with q 1 a real number, and its inverse, the g-deformed exponential function 
exp 9 (x) = [f + (f — q)x] 1 ^ 1 ^ that reduce, respectively, to the ordinary loga- 
rithmic and exponential functions when q = 1. In terms of these functions, Eq. 
(2) and its inverse can be written more economically as 

log Q N(k) = \og a A max - Af-'k, (3) 

and 

N(k) = N^expJ-N^Af-'k}. (4) 

We first comment that Eq. (4) is a generalization of Zipf's law that takes properly 
into account the behavior for low rank k observed in real data where, as one 
would expect, iV max is finite. In Fig. 1 we compare the numbers of occurrences of 
English words in a corpus with N(k) as given by Eq. (4) where the reproduction 
of the small-rank bend displayed by the data before the power-law behavior sets 
in is evident. In the theoretical expression this regime persists up to infinite rank 
k — > oo. Alternatively, we recover from Eq. (4) the power law N(k) ~ fc 1 ^ 1- ") 
in the limit -/V max 3> 1 when a > 1. We note that for ranked listings of data 
N the normalization of their distribution P(N) implies that the maximum rank 
fcmax is equal to the number of data N. Normalization of P(N) = N~ a leads to 
fcmax — A" with both fc max — > oo and J\f — > oo, but J\f~ 1 k generally finite. The 
assumption of a pure power law form for P(N) cannot represent a set with a 
finite number of data. 

Now, in order to arrive at an interesting physical interpretation of Eq. (3) 
we look at the quantities contained in it. We notice that both \og a N max and 
log a N(k) are given by the integrals 

log Q iV max = J N~ a dN and log Q N(k) — J N- a dN, (5) 
l l 

and these in turn can be seen, when a = 1, to conform to the evaluation of 
entropy Si = logAf max or Si — log N(k) where the probability of N equally- 
probable configurations in phase space is P(N) — N~ x . If we now allow for 
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Fig. 1. Rank-order statistics for the occurrence of words (empty circles) in the British 
National Corpus [5]. Eq. (4) with a — 2.09 (smooth curve) is fitted to the data. The 
straight line in the inset is drawn for visualization purposes. 



a > 1 we can retain the same interpretation, 

S a = log a N max and S a = \og a N(k), (6) 

with P(N) = N~ a still viewed as the probability of N equally-probable phase- 
space configurations, and with 7V max and JV(fc) playing the roles of total con- 
figurational numbers or partition functions. Therefore Eq. (3) can be rewritten 

as 

S a = S a -Af~ 1 k, (7) 

and read as the expression of what we refer to as an incomplete Legendre trans- 
form from the Massieu potential S^A/" -1 ), a function of the inverse of the num- 
ber Af, to the entropy S a (k), a function of the rank k. The conjugate variables 
N^ 1 and k could be seen to play the roles, for example, of inverse temperature 
j3 and energy u in the description of a thermal system. As we know the Legendre 
transform is performed in two steps, the first is to add (subtract) the product 
of two conjugate variables from one thermodynamic potential and the second 
is to eliminate the variable in the first potential in favor of the other variable 
to obtain the second potential. The last step involves the derivative of the first 
potential, as the Legendre transform is associated to an extremum value. But 
stopping the procedure at the first step and use of the generalized potential that 
depends on the two conjugate variables is not devoid of use. Familiar examples 
of incomplete Legendre transforms are the Landau free energy (when describing 
a magnet it has a dependence on both magnetization and external field) and 
the free energy density functionals associated to many thermal problems. Eq. 
(4), being the inverse of Eq. (3), states the same relationship but in terms of 
the 'partition functions' N(k) and A r max (A/' _1 ). The absence of an upper bound 
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for the rank k indicates a condition we refer to as the thermodynamic limit in 
our statistical-mechanical interpretation of Eq. (3). To complete the Legendre 
transformation of S a (Af~ 1 ) into S a (k) and eliminate the variable Af^ 1 in favor 
of k, it would be required to optimize S a , i.e. via the use of an 'equation of state' 

fc-^rrlogoAWAT 1 ). (8) 

Wc address this issue in more detail in Section 5. 



4 Analogy with the tangent bifurcation 

Remarkably, there is a strict analogy between the generalized law of Zipf, Eqs. (3) 
and (4), and the nonlinear dynamics for the RG fixed-point map at the tangent 
bifurcation, as originally realized in Ref. [8] . Consequently, these two apparently 
different problems share the same statistical-mechanical interpretation indicated 
in the previous section, and the equivalence offers an alternative to advance our 
analysis, specifically, the characterization of finite size effects for the generalized 
law in terms of the shift of the map out of tangency. 

The analogy can be seen immediately after a brief recall of the RG treatment 
of the tangent bifurcation that mediates the transition between chaotic and 
periodic attractors [6]. The common procedure to study the transition to chaos 
from a trajectory of period n starts with the composition f <<n \x) of a one- 
dimensional map f(x) at such bifurcation, followed by an expansion for the 
neighborhood of one of the n points tangent to the line with unit slope [6] . With 
complete generality one obtains 

x' = f (n) (x) =x + ux z + x > 0, z > 1, (9) 
where x z =sign(a;) \x\ z . The RG fixed-point map is the solution f*(x) of 

/*(/*(*)) = A" 1 /* (A*) (10) 

together with a specific value for A that upon expansion around x = reproduces 
Eq. (9). An exact analytical expression for f*(x) was obtained in Ref. [8] with 
the use of the assumed translation property of an auxiliary variable, y = x x ~ z . 
This property is written as 

x' 1 - 2 = x 1 - 7 ' + (1 - z)u (11) 

or, equivalently, as 

x' = xcx-p^ux*- 1 ). (12) 

It is straightforward to corroborate that x' — f*(x) as given by Eq. (12) satisfies 
Eq. (10) with A = 2 1 ^ z - 1 \ Repeated iteration of Eq. (11) leads to 

x]- z = x\- z + (1 - z)ut (13) 
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or 

log z x t = log z x + ut. (14) 
So that the iteration number or time t dependence of all trajectories is given by 

Xt = x exp z [xq~ ut\ , (15) 

where the xq are the initial positions. The q-deformcd properties of the tangent 
bifurcation are discussed at greater length in Ref. [7]. The parallel between Eqs. 
(14) and (15) with Eqs. (3) and (4), respectively, is plain, and therefore we 
conclude that the dynamical system represented by the fixed-point map f*(x) 
operates in accordance to the same statistical-mechanical property described in 
the previous section for the generalized laws. 

To emphasize that there is a firm analogy, not a casual resemblance, between 
the ranking of data and the sequences of iterates at the tangent bifurcation we 
show that there is a common source behind Eqs. (3) and (14), i.e. the restriction 
of accessibility to phase space already mentioned. This is readily seen by con- 
sidering the replacement, valid for large time r, of the difference x T+ i — x T by 
dx T /dT in Eq. (9), written as x T+ i — x T = u \x T \ z . Integration of the left hand 
side of the resulting differential form 

iff = udT (16) 

between x and x t and the right hand side from to t leads immediately to Eqs. 
(13) or (14). The quantity |£t|~ z in Eq. (16) plays the same role as the power 
law distribution P{N) ~ N~ a . 

We notice that the absence of an upper bound for the rank k in Eqs. (3) 
and (4) is equivalent to the tangency condition in the map. Accordingly, we look 
at the changes in N(k) brought about by shifting the corresponding map from 
tangency (see Fig. 2), i.e., we consider the trajectories x t with initial positions 
x a of the map 

x' = xexp^ux 2 ^ 1 ) + e, < e < 1 (17) 

with the identifications k = t, N^ 1 = —u, N(k) — x t + x* , 7V max = x + x* and 
a = z, where the translation x* ensures that all N(k) > 0. 

In Figs. 3 to 5 we illustrate the capability of this approach to reproduce 
quantitatively real data for ranking of eigenfactors (a measure of the overall 
value) of physics journals [9], industrial production growth rates by country [10], 
and per capita carbon dioxide emissions by country or region [11], respectively. 

In the intermittency route out of chaos it is relevant to determine the du- 
ration of the so-called laminar episodes [6], i.e., the average time spent by the 
trajectories going through the "bottle neck" formed in the region where the 
map is closest to the line of unit slope. Naturally, the duration of the laminar 
episodes diverges at the tangent bifurcation when the Lyapunov exponent for 
separation of trajectories vanishes. Interestingly, it is this property of the non- 
linear dynamics that translates into the finite-size (fc max < oo) properties of the 
occurrence-rank function N(k), that we have obtained without finding out the 
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Fig. 2. The map in Eq.(17) with a trajectory. The inset shows the time dependence of 
the trajectory. 

details of the departure of the basic distribution P(N) from the pure power- 
law N~ a . One more important result that follows from the analogy between 
nonlinear dynamics and the rank law is that the most common value for the 
degree of nonlinearity at tangency is z — 2, obtained when the map is analytic 
at x = with nonzero second derivative, and this implies a = 2, close to the 
values observed for most sets of real data. 

5 Universality and uniqueness of data ranking 

The steepest-descent approximation is central to statistical mechanics (and in a 
more general context to large deviation theory [13]). This property facilitates the 
evaluation in the thermodynamic limit of a partition function for one particular 
ensemble in terms of the partition function of another. Thermodynamically, this 
approximation relates to the Legendre transformation between the corresponding 
free energies or Massieu potentials where one variable is eliminated in favor of its 
conjugate [14]. As recalled above, the procedure consists of two steps, summation 
(or subtraction) of the product of the conjugate variables to (or from) the first 
potential to define the second, followed by use of the derivative of the first 
potential, or equation of state, to remove the undesired variable. This, of course, 
corresponds to the optimization involved in the steepest-descent method. For 
illustrative purposes we will assume here that the steepest-descent shortcut that 
underlies the second step in the Legendre transformation is also meaningful for 
a > 1. 

In order to carry out the second step in the Legendre transform stated by 
Eq. (7) we need an explicit form for the function A r max (A/' _1 ). It is evident that 
the form of this function is not unique and is determined by the particular set 
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Fig. 4. Rank-order statistics for industrial production growth rates (empty circles) 
from Ref. [10]. Eq. (17) with the identifications provided in the text when a = 2.13 
and e = —0.058 (smooth curve) is fitted to the data. 
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Fig. 5. Rank-order statistics for per capita carbon dioxide emissions (empty circles) 
from Rcf.[ll]. Eq. (17) with the identifications provided in the text whena = 2.06 and 
e — —0.0055 (smooth curve) is fitted to the data. 

of data Af. For illustrative purposes we consider a finite set of data N extracted 
from P(N) = N~ a although a pure power law is not the correct distribution 
in this case. However, if the equivalent map Eq. (17) is very close to tangency 
£«1 and the data for the maximum rank, iV m i n = N(k max ), is chosen such that 
its image in the map is to the left and near to its bottleneck mid-point, then 
P(N) is closely approximated by the power lawiV~ Q . Under this approximation 
normalization of P(N) only yields k max ~ Af. Suppose the available data, or the 
choice of the data collector, fixes the specific value of k max and the lowest and 
upper limits in Eq. (2) to be N m i D and N max , respectively. Therefore we have 

N max 

*W=AA J N~ a dN=^-[N^-N^], (18) 

a 7^ 1, or 

iV max = AU„ expJiV^AT i/w]. (19) 

For example, a set of data about population of cities may be represented by 
fcmax = 50 (fifty representative city sizes), -/V min = 1 (one city with the largest 
population), and N max = 100 (one hundred cities with the smallest population 
considered). Eq. (19) is the required expression for iV max (A/" _1 ) to be used in 
the 'steepest-descent condition' or 'equation of state' Eq. (8). The result follows 
immediately, it is k = fc max . 

As in ordinary thermodynamics, we observe that the universality of the laws 
described by Eqs. (3) and (4) is due to the general form of the incomplete 
Legendre transformation, while the specific forms adopted by the potentials S a = 
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log a -/V max (.A/" _1 ) and S a = log Q N(k) are particular to the system or situation 
considered. 

6 Summary and discussion 

We have suggested here a novel thermodynamic, or statistical-mechanical, in- 
terpretation or understanding of the generalized laws of Benford and Zipf. The 
expressions for these laws, Eqs. (1) and (2) (or alternatively (4)) were derived 
in Ref. [4] under the basic assumption that the data sets obeyed by these laws 
are statistically well reproduced when extracted from a power law distribution 
P(N) <~ N~ a . We remark here that the deviation from unity of the exponent 
a implies a restricted access to the phase space for the data configurations that 
when enumerated produce the numbers N. The restriction involves an accessible 
subset of this space with a scale invariant property, i.e., a fractal set, as implied 
by the power law N~ a . This viewpoint becomes evident when P(N) is seen to 
represent the probability distribution of N equally probable configurations in 
the phase space for the data, and, consequently, suggests the definition of the 
generalized entropies in Eq. (6). It is important to clarify that the statistical- 
mechanical structure considered here and obtained from the usual via a scalar 
deformation parameter (represented by the power a) docs not conform to that 
known as nonextensive statistics [15] [16]. Even though we define entropies or 
Massieu potentials with the use of the g-logarithmic function and make use of its 
inverse, the g-exponential, we do not require or implicate the optimization of any 
of these quantities via the use of the constraints employed in the nonextensive 
formalism or involve the use of the so-called escort distributions [16]. 

The ranking of real data habitually shows deviations from the Zipf 's power- 
law regime both for small and large rank that can be clearly observed in semi-log 
plots. As we have shown in Fig. 1 the generalized Zipf 's law given by Eq. (4) is 
capable of reproducing accurately the low rank deviation but not that for large 
rank as the power-law regime in this equation extends to k — > oo. An upper 
bound for k suggests finite-size effects inherent in real data. We have captured 
the nature of the upper bound for k by first demonstrating a precise analogy 
between the expression for the ranking laws, Eqs. (3) and (4), and those for the 
dynamics at the transition to chaos via intermittency (the tangent bifurcation) 
in nonlinear maps of low dimensions. The finite-size effects in the ranking of data 
are seen to correspond to the shift off tangency in the map, so that the position of 
the upper bound for the rank k is given by the duration of the laminar episodes 
of chaotic trajectories near the transition to regular behavior. Interestingly, the 
statistical-mechanical interpretation put forward for the generalized law of Zipf 
extends over to the critical dynamics of the transition to chaos via intermittency. 
While, on the practical side, data for the ranking of data for all k is reproduced 
quantitatively by our formalism, as illustrated, respectively, in Figs. 3 to 5 for 
three specific examples: cigenfactor of physics journals [9], industrial production 
rates [10], and carbon emissions [11]. In agreement with empirical determinations 
the analogy implies that the most general value for the index a is a = 2. 
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As it is generally well-known, a statistical-mechanical structure (shared by 
large deviation theory [13]) is built around the steepest-descent approximation 
and is expressed as the Legendre transform property that links different thermo- 
dynamic potentials. It involves an optimization condition or equation of state 
that relates conjugate variables. Only for illustrative purposes we have assumed 
that this structure extends to the deformed version (with one scalar parameter) 
we have considered here. In order to replicate the circumstances normally en- 
countered in thermodynamics we have presented as an example the particular 
form taken by the function iV max (A/'~ 1 ) when the data in hand is bounded by the 
numbers N m { n and iV max and fixes the largest rank fc max . Then, the equation of 
state was determined and the variable A/" -1 eliminated in favor of k, to obtain the 
'equilibrium' value for N(k). This exercise suggests that the universality of the 
laws is due to the general form of the incomplete Legendre transformation, while 
the expressions for the initial and transformed potentials are specific to the de- 
sign of the data sample under consideration. The thermodynamic interpretation 
we have put forward may explain the ever presence of these phenomenological 
laws in a wide range of observations including very dissimilar situations. Finally, 
we comment that our arguments also apply to the topic of scale-free networks 
[17]. Since the degree distribution p(k), the distribution for the number k of links 
that connect one node to other nodes, describes essentially the ranking of nodes 
according to the number of links they possess, we can treat the data sets from 
where this distribution is phenomenologically obtained similarly to the data sets 
leading to Zipf's law. Interestingly, for random link networks p(k) decays ex- 
ponentially (a = 1), but for scale- free networks it is approximately power law 
(a>l). ' 
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